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ABSTRACT 

One model for the origin of typical galactic star clusters such as the Orion Nebula Cluster (ONC) is 
that they form via the rapid, efficient collapse of a bound gas clump within a larger, gravitationally- 
unbound giant molecular cloud. However, simulations in support of this scenario have thus far have 
not included the radiation feedback produced by the stars; radiative simulations have been limited 
to significantly smaller or lower density regions. Here we use the ORION adaptive mesh refinement 
code to conduct the first ever radiation-hydrodynamic simulations of the global collapse scenario for 
the formation of an ONC-like cluster. We show that radiative feedback has a dramatic effect on the 
evolution: once the first ~ 10 — 20% of the gas mass is incorporated into stars, their radiative feedback 
raises the gas temperature high enough to suppress any further fragmentation. However, gas continues 
to accrete onto existing stars, and, as a result, the stellar mass distribution becomes increasingly 
top-heavy, eventually rendering it incompatible with the observed IMF. Systematic variation in the 
location of the IMF peak as star formation proceeds is incompatible with the observed invariance of 
the IMF between star clusters, unless some unknown mechanism synchronizes the IMFs in different 
clusters by ensuring that star formation is always truncated when the IMF peak reaches a particular 
value. We therefore conclude that the global collapse scenario, at least in its simplest form, is not 
compatible with the observed stellar IMF. We speculate that processes that slow down star formation, 
and thus reduce the accretion luminosity, may be able to resolve the problem. 

Subject headings: ISM: clouds — radiative transfer — stars: formation — stars: luminosity function, 
mass function — turbulence 



1. INTRODUCTION 



The origin of the stellar initial mass function (IMF) 
is one of the outstanding problems in the modern the- 
ory of star formation. While there have been numerous 
analytic and numerical studies purporting to explain its 



origin (e.g. see the review by McKee & Ostriker 2007 
and references therein) , much of this work has been ham- 
pered by the limited number of physical processes that 
are included in models of how gas fragments. In par- 
ticular, while both simulations and analytic work reveal 
that how gas fragments into stars is extremely sensitive 
to how the temperature of t he gas varies with its density 
( Larson|[2005 Jappsen et al.| |2005), it has been common 
until very recently to approximate this relationship with 
a simple equation of state (e.g. Bate & Bonnell 2005; 
Bonneh et al. 120061 [Offner et al..2008b..Hennebelle et al.. 
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201ip . Since the characteristic masses of the stars formed 
m a collapse are largely determined by the temperature- 
density relationship, predictions about the location of the 
IMF peak in these simulations are only as good as their 
adopted equations of state. 

Given this realization, attention in recent years 
has shifted to models that attempt to determine the 
temperature-density relationship from first principles, or 
to include a self-consistent treatment of the thermal evo- 
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lution of the gas in numerical simulations. In the former 
category, much work has focused on the effects of imper- 
fect coupling between gas and dus t grains on gas thermo- 
dynamics. For example, [Larson] ([2005 ) and Elmeg reen| 
|et al.| ( |2008|) both argue that the characteristic stellar 
mass is set by the Jeans mass at the density and tem- 
perature where dust grains and gas become thermally 
coupled due to collisions. According to these models, at 
low densities where grain-gas coupling is poor, the gas 
is slightly sub-isothermal, while at higher densities it is 
slightly super-isothermal, and this effect favors fragmen- 
tation near the coupling density. 

However, this argument faces a major difficulty in ex- 
plaining the IMF in the dense, cluster- forming regions 
where much Galactic star formation appears to take 
place. The density at which grains and gas bec ome well- 
coupl ed is ~ 10^ — 10^ H2 molecules cm~^ ([Goldsmith 



2001), roughly independent of the m etallicit y and of aiii^ 
bient radiation field intensity ([Krum holz & McKee|2008 



Elmegreen et al. 2008 Krumholz et al.n201l| ). In com 



parison, observations now show that the typical site of 
star cluster formation has a mass of ^ 10^ — 10^ M^, and 



a radius of ^ 0.3 — 0.5 pc (e.g. see | Shirley et al. 



2003 



Faundez et al. 120041 |Fontani et al. 2005 , or the su mmary 



plot combining these data sets in , Fall et al.|2010 ), giving 
a mean density ~ 10^ cm~^. Similarly, the present-day 
Orion Nebula Cluster has a mass of 2400 Mq within 
a half mass radius of 0.8 pc, corresponding to 2 x 10^ 
cm~^, and within the ^ 0.2 pc core the mean density 
reaches 4 x 10^ cm~^ (Hillenbrand & Hartmann 1998). 



Since the star formation etiiciency was certainly less than 
unity, and the cluster has likely expanded some since the 
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gas was expelled rtKroupa et al. 2001 Tan et al. 2006L 
the density at which most of the stars formed must have 
been higher by at least a factor of a few. Thus the typical 
site of star cluster formation in the Galaxy, of which the 
ONC is an example, is in the regime where essentially 
all the mass is at densities where grain-gas coupling is 
very strong. It is therefore hard to see how grain-gas 
coupling could be relevant for determining how this gas 
fragments. This argument can be made even stronger by 
noting that globular clusters with mean densities ~ 10^ 
cm~^ in their centers, 2 — 3 orders of magnitude above 
the grain-gas decoupling density, also appear t o have the 



same IMF peak as the Galactic field ( Marchesini et al. 
|2009). 

A second class of models for the temperature-density 
relationship focuses on the interaction of gas with the ra- 
diation produced by stars in the star formation process. 
In these models one assumes good grain-gas coupling, 
as is appropriate at the high densities where most stars 
form. The gas temperature and its relationship with the 
density is then determined primarily by the light pro- 
duced by stars in the process of formation. Conceptu- 
ally, the idea is that the luminosity from an accreting 
star warms the gas in its immediate vicinity, inhibiting 
the ability of that gas to fragment, and that this pro- 
cess determines charac teristic stellar mas ses. Analyti- 
cally, |Krumholz| (j2006) and K rumholz fc McKee| ( |2008 ) 
have argued that this process explains how massive stars 
are able to form under certain circumstances, while Bate 
|2009) argues that it can explain the characteristic peak 
of the IMF. 

However, numerical studies of the second class of mod- 



els have thus far be en limited i n various wa ys. [Krumholz 
|et al.| ( |2007a[ [2010) and Myers et alj ( |2011l ) 
ulations including; stellar feedback and raai 



conduct sim- 
liative transfer 
(including re-radiation by dust grains, which is the crit- 
ical process in determining the gas temperature), but 
focus on single massive cores that do not (and are not 
expected to) form a full IMF. Commercon et al. ( [2010 ) 
report similar simul ations focusing on si ngle low-mass 



report simnar simul ations locusmg on si ngle low-mass 
cores. [Batel (|2009|, |Offne^_e^ aPjgOOQl) , [Price & Bate] 
(I2009|), ana nFeters etal.| ( |2010[ |201l|) simulate the for- 



mation of star clusters, but consider only low-density re- 
gions similar to those found in nearby clouds like Taurus, 
rather than conditions typical of Galactic star formation 
sitesj^ As we will see, this makes a large difference in 
the outcome, because under low-density conditions the 
regions of heating around each star are non-overlapping, 
while in denser condit ions they are no t. Moreover, of 



these simulations, only [Off ner et al. and Peters et al. 



in- 
clude stellar luminosity, so the amount of heating in the 
other two simulations is underestimated. 

Other simulations of star cluster formation do not in- 
clude radiative transfer at all, a nd instead ap proximate 
it in various ways. For example, Smith et al. (2009) and 



Urban et al. (2010) study the fragmentation of dense gas 



^ Although Peters et al.] ( |2010[[20TT] > study regions with enough 
mass to form massive stars, the coluixin densities of the regions 
they simulate are ~ 0.01 g cm"-^, rather than the ~ 1 g cm~^ 
typical of Galactic star-forming sites. Their simulated clouds are 
therefore optically thin even in the near-IR, rendering radiative 
effects fairly unimportant . In this way t heir work is closer to that 
of |Bate[|Offner et al.( and [Price &: Bate] than to the simulations we 
present here! 



clouds similar to typical star-forming regions, but they 
determine the gas temperature around each star via a 
rough fitting formula based on static radiative transfer 
calculations. This approximation may be reasonable as 
long as the heating at a given point is dominated by a 
single star, but it almost certainly fails once the regions 
of heating around stars begin to overlap, as occurs in 
dense regions. 

In summary, to date there have been no simulations ca- 
pable of studying how the peak of the IMF under typical 
Galactic conditions, including the all-important effects of 
stellar feedback and re-radiation by dust grains. The goal 
of this paper, the first in a series, is to remedy that lack. 
We use the ORION adaptive mesh refinement radiation- 
hydrodynamics code to simulate a typical galactic star- 
forming clump including stellar feedback and radiative 
transfer. As this is a first attack on the problem, we 
choose the simplest possible scenario. We do not include 
magnetic fields or any form of feedback other than radi- 
ation, and we allow the initial turbulence in the cloud to 
decay freely, leading to a rapid global collapse. Our sim- 
ulation therefore represents a minimalist scenario for the 
formation of a star cluster such as the ONC similar to 



that proposed by, for example Bonneh et al.| ( (2003[ ). Pre- 
vious authors who have studied such conditions report 
that they produce stellar mass distributions consistent 
with the observed IMF at all times in the simulation, 
but it is clear in retrospect that this result simply re- 
flects the imposed equation of state. Our work therefore 
revisits the critical question of whether such a scenario 
is capable of reproducing the observed IMF. 

The remainder of this paper proceeds as follows. In 
Section [2] we describe our numerical method and simu- 
lation setup. In Section [3| we report the results of our 
simulations. In Section [ijwe discuss the implications of 
our findings, and present simple analytic models to aid in 
understanding them. Finally, we summarize in Section 

El 

2. SIMULATION DESCRIPTION 
2.1. Simulation Initial Conditions 

We conduct two simulations that are identical in every 
respect except that they have different maximum AMR 
levels, meaning that the peak resolution is different in the 
two runs. We refer to these as the low-resolution (LR) 
and high-resolution (HR) simulations. The two simula- 
tions enable us to determine to what extent our results 
are converged, although we caution that the two simu- 
lations differ only in their peak resolution, which is de- 
ployed near stars and in regions of high density or large 
radiation gradients (see Section 2.4). Thus we have not 
tested the sensitivity of our results to variations in the 
resolution used in low density regions far from stars. We 
also carry out a third simulation with identical initial 
conditions at the same resolution as run HR, but with 
an isothermal equation of state, i.e. with the radiative 
transfer module in our code disabled. We refer to this as 
run ISO. This simulation enables us to determine what 
effects in our simulation are due to radiative transfer. 

We summarize the key parameters of the runs in Table 
[l] The initial conditions for both consist of a Mc = 1000 
Mq spherical gas cloud with a mean surface density He = 
1 g cm~^, corresponding to a mean volume density of 
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TABLE 1 
Simulations 



Name 


RT? 


Mc (Mo) 


Ec (g cm-2) 


p (g cm 


') 


(kyr) 


^box (pc) 




L 


AxL (AU) 






LR 


Yes 


1000 


1.0 


9.4 X 10- 


-19 


68.6 


1.9 


256 


4 


98 


0.94 


0.51 


HR 


Yes 


1000 


1.0 


9.4 X 10- 


-19 


68.6 


1.9 


256 


5 


49 


0.94 


0.52 


ISO 


No 


1000 


1.0 


9.4 X 10- 


-19 


68.6 


1.9 


256 


5 


49 


0.94 


0.65 



Note. — Col. 2: radiative transfer included? Col. 3: cloud mass. Col. 4: cloud surface density. Col. 5: mean volume 
density. Col. 6: mean-density free-fall time. Col. 7: linear size of computational domain. Col. 8: number of grid cells per linear 
dimension on the coarsest AMR level. Col. 9: maximum AMR level. Col. 10: linear cell size on the maximum AMR level. 
Col. 11: time relative to free-fall time to which simulation is evolved. Col. 12: total mass of stars at the final evolution time, 
normalized to the initial cloud mass. 



9.4 X 10"^^ g cm-^, or 2.4 x 10^ H2 mo lecules cm~ ^. The 
corresponding cloud radius is Rc = \/Mc/ (ttSc) =0.26 
pc, and we place the cloud in a cubical computational 
domain of size Lbox = 1-9 pc, roughly four times the 
cloud diameter. We have chosen this mass and surface 
density because they are typical of regions of clus tered 
star fo rmation in th e Galaxy (e.g. see Shirley et al.|2QQ3 , 
'Faundez et al.]|2QQ4[ [Fontani et al. ..20051 Q-nd a sunimary 
of the data in Figure 1 of |Fan et al.||2QTQ| ). They are 



also roughly the estimated parameters of t he progenitor 
of the core of the Orion Nebula Cluster (e.g. |Kroupa et alT 
2001 Tan et al. 2006). It is worth noting that our initial 
conditions are significantly denser than has been used for 
some pre vious simulations of m assive star formation. For 
example, Bonnell et al. (2003) use initial mean volume 
and column densities of 1.3 x 10~^^ g cm~^ (3.3 x 10^ 



cm and 0.26 g cm , respectively;^ 

use3.9xl0"^^ gcm-^ (1.0x10^ cm"^) and 0.026 g cnT^. 
However, our parameter choices are much closer to what 
is actually observed in regions of massive star formation. 
For example, in thei r survey of 1 46 Southern massive 
star- forming regions, Faundez et alj ([2004j) find a typical 
mass and radius of 5000 Mq and 0.4 pc, corresponding 
to a volume density of 1.2 x 10~^^ g cm~^ (3.1 x 10^ 
cm~^) and a column density of 2.1 g cm~^, similar to 
what we use. 

Our initial cloud has a density structure described by 



^JGMc/2Rc = 2.9 km s ^. The corresponding virial 
parameter is a = ba^GMc/ Rc = 5, so that the turbu- 
lent kinetic energy is larger than the potential energy at 
time zero. However, we do not include any feedback pro- 
cesses (e.g. winds or H ll regions) capable of driving the 
turbulence, nor do we have other potential driving mech- 
anisms, such as a turbulent cascade from larger scales or 
ongoing infall. As a result, the turbulence undergoes a 
rapid decay, which quickly renders the cloud gravitation- 
ally bound. 

Throughout the computational domain, we initialize 
the radiation energy density to that of a blackbody with 
a temperature = 10 K. Thus we have E = aT^ = 
7.56 X 10~^^ erg cm~^. Similarly, we initialize the gas 
temperature within the cloud (r < Rc) to Tg = 10 K. 

Outside the cloud (r > Rc)^ we set the temperature to 

Peters et al.|(|2010|) Tg = 1000 K. Since the density outside the cloud is 1/100 



Pc, 

Pc{2r/Rc)-^-^ 
2-i-V,/100, 



r < Rc/2 
Rc/2 <r <Rc 
r>Rc 



that of the density at the cloud edge, this ensures thermal 
pressure balance across the cloud boundary. We also 
set the Planck and Rosseland opacities of the material 
with Tg > 500 K and p < 2~^-^pc/50 to zero, to ensure 
that the host ambient medium does not interact with the 
radiation field, and is not able to cool. 

2.2. Evolution Equations 

The simulations we present in this paper use the 
ORION adaptive mesh refinement code. The numeri- 
cal method is nearly identical to that in our previous 



papers (e.g.JKr umholz et al."2007al |2009| |2010[ [Myers 
iCuhni 



(^l^ et al.|20lT Cunningham et al._2011) . Here we only sum 



where r is the distance from the cloud center and pc = 
6^c/[{2^'^ - l)^c] = 1.6 X 10-^^ g cm-3 is the core den- 
sity. Thus our density profile consists of a constant den- 
sity in the inner half of the radius, coupled with a pow- 
erlaw falloff in the outer half of the radius. Outside this 
cloud we place a low density ambient medium with a den- 
sity that is 100 times smaller than the cloud edge density. 
We choose this density structure because observations in- 
dicate the presence of a roughly r~^-^ densi ty gradient on 



large-scal es in star -forming clu mps (e.g. iCaselli fc Myers 
1995; Be uther et al . 2002, 20051 [20061 [Mueller et al.[M7^ 
Sridhara n et al.[[2Q05, ). By choosing a fiat inner density 
profile, however, we minimize tidal forces in the cloud 
core, thereby ensuring maximum opportunity for frag- 
mentation. 

We initialize the cloud velocity with a Gaussian- 
random velocity field with a power spectrum P{k) ex 
k~'^ and a one-dimensional velocity dispersion ctc = 



marize the physics, and we reter rea ders to the numerical 
method papers referenced in Section [23] for a full descrip- 
tion of ORION's workings. ORION works by solving the 
equations of compressible gas dynamics including self- 
gravity, radiative transfer, and radiating star particles, 
all on an adaptive grid. In our computational domain, 
we describe every cell with a vector of conserved quan- 
tities (p, pv, pe, where p is the density, pv is the mo- 
mentum density, pe is the total internal plus kinetic gas 
energy density, and E is the radiation energy density in 
the rest frame of the computational grid. In addition 
to gas quantities, we also track an arbitrary number of 
point mass star particles, each of which is described by 
a position cc^, a momentum p^, and an instantaneous lu- 
minosity L^, where the subscript i refers to the particle 
number. 

Given this description of the problem, the full set of 
evolution equations is 



d 



■{pv)-Y,M^W{x-X,) 



(2) 
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dt 



(PV): 



-V • {pvv) - VP - pV(f> - WE 



(3) 



dt 



ipe): 



-V • [{pe + P)v] - pv-V(f>- KoppiinB - cE) 



Xl2'^-l]vVE 



dt 
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cX 



P 



rrir 



(4) 
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yE + Kopp{4:nB - cE) 



l]vVE-V-{ ^-J^vE 



dt 



Mi = M, 



dt 
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di^ 



Xi - 



Pi. 

Mi 



(-^) A{T,) + J2LiW{x-Xi) (5) 

(6) 
(7) 
(8) 



-MiV0 + Pi 



-AttG 



p + J^^i^i^-^) 



(9) 



Equations ([3|, and Q represent the conservation 
laws for gas mass, momentum, and energy, including 
terms describing the exchange of these quantities with 
star particles and with the radiation field. Equation ^ is 
the corresponding conservation of energy equation for the 
radiation field. Similarly, equations (|6|, ([7|), and ([8| are 
the equations of mass and momentum conservation, and 
the equation of motion, for the point particles. Finally, 
equation ^ is the Poisson equation for the gravitational 
potential Note that we compute the gas-radiation 
exchange terms using the mixed frame formulation ( |Mi-| 
halas fc Klein 1982; Krumholz et al. 2007b), allowing us 
to write them in a form that is manifestly and exactly 
energy- conserving . 

In these equations, the terms M^, p^, and Ei represent 
the rate at which mass, momentum, and energy accrete 
from the gas onto the ith star, and Li represents the 
luminosity of that star. We de scribe how we compute 
these quantities in Section 2.3 The quantities P and 
Tg are the pressure and gas temperature, respectively. 
These are related by the equation of state 



P = 



(10) 



where /i = 2.33 is the mean molecular weight for molecu- 
lar gas of Solar composition and 7 is the ratio of specific 
heats. We adopt 7 = 5/3, appropriate for gas too cool 
for hydrogen to be rotationally excited, but this choice 
is essentially irrelevant because Tg is set almost purely 
by radiative effects. The quantities hcqp and a^or are the 
specific Planck- and Rosseland-mean opacities in the rest 
frame of the gas, B = caRT^/An is the Planck function. 



and A is the fiux limiter, given by 



R-- 



\VE\ 



KORpE 

R2=\ + \^R^. 



(11) 

(12) 
(13) 

We compute the opacities as a function of the gas den- 
sity and temperature using th e iron normal, c omposite 
aggregates dust model of Semenov et al. (2003). 

Finally, A(T^) represents the line cooling coefficient. 
We include this because the turbulence can be strong 
enough in our cluster so that, at isolated points, gas 
shock- heats to temperatures above a few thousand K. 
This exceeds the dust sublimation temperature, so the 
dust opacity becomes nearly zero in this gas. Instead, 
the gas in this temperature regime is cooled by line emis- 
sion, which we cannot easily describe with a simple con- 
tinuum opacity. In this gas, we transfer energy from 
the gas thermal reservoir to the radiation field at a rate 
{p/rripY K{Tg)^ where rrip is the proton mass, and th e 



fun ction K{Ta) is taken from Cunningham et al. (2006). 



See [Cunningham et al. (2011 ) tor more details ot our line 
cooling approach. 

An important subtlety in our evolution equations, 
which is worth noting, is that we do not differentiate 
between gas and dust grain temperatures. At low densi- 
ties, gas-grain coupling can be imperfect, and it can be 
important to calculate the two temperatures separately, 
and to s imulate the thermal exchange between dust and 
gas (e.g. Urban et al.|2009[ ). However, grains and gas be- 
come very tightly coupled in temperature once the den- 



sity exceeds n 



10^ 



10 cm . For comparison, the 



mean density in our initial clouds is n = p/ja = 4.0 x 10^ 
cm~^. Thus our entire computation is in the strong cou- 
pling regime, and there is no need to treat dust and gas 
temperatures separately. 

For simulation ISO, our isothermal run, we modify 
these equations as follows. First, we omit equation ^ 
entirely. Second we set to zero all terms proportional to 
E or A{Tg) in equations (3| and (|4|. Third, instead of 
7 = 5/3, we use 7 = 1. 000 iT This corresponds to neglect- 
ing the effects of radiative transfer, and simply keeping 
the gas temperature almost completely fixed to its initial 
value. 

2.3. Numerical Method 

The ORION code solves equations ([2| - ([9| in a series 
of operator-split steps. In each time step, we first inte- 
grate the hydrodynamic equations ([2| - (|4|, excluding 
the terms describing stars and the radiation field. This 
update uses a conservative Godunov scheme with an ap- 
proximate Riemann solver, and is second-order accurate 



in time and space (Truelove et al. 1998 Klein 1999) 



Next we solve the ro isson equation j9| using a multi- 



grid ite ration scheme (Truelove et al. T598[ |Klein|p!999 
Fisher, 2002). Third, in the runs where we include radia- 



tion, we update the radiation energy equation (|5| and the 
radiation terms in the hydrodynamic equation s ([2| - (|9|. 
This update uses the Krumholz et al. ( 2007b[ ) conserva- 
tive update scheme, m which we handle the dominant 
terms implicitly and the non-dominant terms explicitly. 
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The update f or the imphcit terms uses the Shestakov fc| 
Offner (2008) pseudo-transient continuation scheme. Fi- 
nahy, we update the stehar quantities, equations (|6| - 
([s]), and update gas quantities for the gas-star exchange 
terms on the right hand sides of equations Q - (|9|. We 
determine the accretion rates of mass, momentum, and 
energy onto each star by fitting the flow within a ra- 
dius of four finest-level cells of each star particle to a 
Bondi-Hoyle flow , follow ing the procedure described by 
Krumholz et al. (2004). We update the luminosity Li 
ot each star using the protostellar evolution model de- 
scribed in the appendices of Offner et al. (2009). 
Each of these update modules operates withi n the over- 



all adaptive mesh framework of ORION (Berger & Oliger 



1984 
schema 



Krumholz et al. 



2007b) and setting 



Berger fc CoIlela][l989l |Bell et al | |1994| ). in this 
e, we discretize the computational domain onto a 
series of levels / = 0,1,2,...,L. The coarsest level, level 
0, has cells of linear size Axq, and covers the entire com- 
putational domain. All subsequent levels, with cells of 
size Axi = Axo/2^, cover subregions of the computa- 
tional domain. Each level consists of a union of rect- 
angular grids of cells, and grids on different levels are 
nested such that every level / grid with / > is fully 
contained within one ore more level / — 1 grids. To ad- 
vance a level / in time, we flrst advance all the grids on 
that level through a time step At^, then advance grids 
on level / + 1 by two timesteps of size At/+i = Ati/2. 
After the two level / + 1 advances, we synchronize the 
boundaries between levels / and / + 1 to ensure exact con- 
servation of mass, momentum, and energy across level 
boundaries. The entire update procedure is recursive, 
so a single advance on level / + 1 entails two advances 
of size = At^+i/2 on level / + 2, and so forth to 

the flnest level present. The coarse level timestep Ato 
is set by computing the Courant condition on each level 
(including a contri bution to the signal sp eed from ra- 
diation pressure - 
Ato = min(2^At/). 

2.4. Boundary, Refinement, and Star Particle 
Conditions 

At the edge of the computational domain, we use 
reflecting boundary conditions for the hydrodynamics. 
However, this choice is irrelevant to the evolution, be- 
cause our computational domain is large enough to en- 
sure that no material from the cloud ever approaches 
it. For the gravity, we adopt Dirichlet boundary condi- 
tions, with the potential at the computational domain 
boundary set equal to a multipole expansion of the po- 
tential due to the matter in the domain interior, includ- 
ing terms up to the quadrupole. Finally, for radiation 
we adopt Marshak boundary conditions, with the flux 
into the computational domain set equal to the flux of 
an isotropic 10 K blackbody: Fin = caT^/A = 0.57 erg 
cm~^ s~^. The boundary condition is equivalent to al- 
lowing any radiation generated within the computational 
domain to escape freely, but also to bathing the compu- 
tational domain in a 10 K blackbody radiation fleld. 

In order to determine when AMR levels are added 
or removed, we must also specify reflnement conditions. 
The conditions we use in our simulations are as follows. 
First, we reflne any cell with a density greater than half 
the edge density of the initial cloud to at least level 1. 



This ensures that our initial cloud is well resolved. Sec- 
ond, we reflne any cell on level / that is within a distance 
16Ax/ or less from any star particle. This ensures that 
the region around each star is resolved by at least 32 cells 
on all levels by the finest one. Third, we refine any cell 
where the density exceeds the local Jeans density ( |Tru-| 
elove et al.|[l997D , 



p> pj = J^ 



GAxf' 



(14) 



where Cg = ^Jk^T^J]l is the sound speed. We use a 
Jeans number J = 1/8. Finally, we refine any cell where 
the local radiation energy gradient satisfied the condi- 
tion \AE\lE > 0.15/Ax/. This ensures that gradients of 
radiation energy density are always well-resolved. If any 
of these conditions are met, we refine that point in the 
computational domain to a higher AMR level, up to the 
maximum level L for that simulation (see Table T]). 

Finally, we create a new star in any cell on the maxi- 
mum level L that violates the Jeans condition, equation 
(14), using a Jeans number J = 1/4. In contrast to pre- 



vious runs, where we merged star particles together if 
they approached closer than a certain limit, here we do 
not allow any star particle that has a mass greater than 
0.05 Mq to be destroyed by merging. Our motivation 
for choosing this mass limit is that it is roughly the mass 
at which sec ond collapse to stellar densities occurs (e.g. 
Masunaga et al . 1998; Masunaga & Inutsuka 2000|). Ob- 
jects of lower mass remain extended gas balls with phys- 
ical sizes of a few AU, and thus are much more likely to 
merge than the much smaller, more compact protostars 
they become once they complete their collapse. Com- 
plete suppression of mergers for more massive objects 
is probably an extreme assumption, but as we will see 
in discussion our results, allowing mergers would only 
strengthen our conclusions, by moving the stellar mass 
distribution to higher values. 

3. SIMULATION RESULTS 

For convenience, throughout this section we will re- 
port our results in terms of mean-density free-fall times, 
where the mean density is p = 3Mc/(47ri?^) = 9.4 x 10"^^ 
g cm~^ and the corresponding free-fall time is tff = 
v/37r/32Gp = 68.6 kyr. The free-fah time in the high- 
density initial core is ~ 30% shorter, tf^^c = 52.3 kyr. In 
reporting stellar quantities, we only count as stars those 
star particles with masses above 0.05 M©, the mass at 
which second collapse to stellar dimensions occurs. How- 
ever, this has little effect on our results, since objects be- 
low this mass never constitute more than a tiny fraction 
of the total mass in star particles. 

We ran these simulations on a combination of the su- 
percomputers Pleiades at the NASA Advanced Super- 
computing facility and Ranger at the Texas Advanced 
Computing Center. Runs LR, HR, and ISO required 
roughly 200,000, 850,000, and 60,000 CPU hours, respec- 
tively, and ran on between 256 and 960 CPUs (32 to 120 
nodes), with the number of CPUs used increasing as a 
run progressed and the number of high resolution grids 
increased. 

3.1. Large-Scale Evolution 
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Fig. 1. — Column density (left) and density- weighted mean tem- 
perature (right) in run LR. The left and right columns show the 
state of the simulation at the same time, with the time running 
from t/tff = at the top to t/% = 0.9 at the bottom, at intervals 
of 0.15. The ratio of stellar mass to initial cloud mass, M*,tot/^c, 
is also indicated in each panel. In the column density plot, white 
circles show the locations of star particles, with the size of the circle 
indicating the mass of the star. In the right column, the tempera- 
ture we show is the radiation temperature, defined by E = a^T^. 
We show this quantity because the radiation and gas temperatures 
are nearly identical everywhere except in the low-density, zero- 
opacity ambient medium, where the radiation temperature is far 
lower than the gas temperature. By using the radiation rather 
than the gas temperature, we ensure that our projected temper- 
atures are not artificially enhanced by contributions from the hot 
ambient medium. 




-0.3 -0.2 -0.1 -0.0 0.1 0.2 -0.3 -0.2 -0.1 -0.0 0.1 0.2 0.3 
X [pc] X [pc] 



Fig. 2. — Same as Figure [T] but for run HR. 

Figures [l] [2j and [3] show the large-scale evolution of 
the cloud as it collapses in our simulations. As the plots 
make clear, the overall distribution of the gas and stellar 
mass, and the gas temperature structure, is very similar 
in all runs. As we will see in more quantitative detail 
later, the evolution of the two radiative runs is very sim- 
ilar in almost every respect, so that we may have confi- 
dence that the behavior we are seeing is physical and not 
a result of resolution effects. Even at late times, the only 
noticeable difference is the exact positions of individual 
stars on the periphery of the cloud. These differ primar- 
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Fig. 3. — Same as Figure p] but for run ISO. Since this run 
is isothermal, we show only column density, not density-weighted 
mean temperature. 



ily because the N-body interactions that occur late in the 
simulation are chaotic. They can therefore be changed 
significantly because the amount of gravitational soften- 
ing in the gas-particle and particle-particle interactions 
is resolution-dependent. 

In both radiative runs, we see that, for the first 
(0.3 — 0.4)tff, the initial velocity perturbations we have 
injected are developing and creating structure, but that 
no stars have yet formed. The gas temperature remains 
locked at 10 K, the value imposed by the radiation field. 
Around t/tff = 0.45, the first stars start to appear at the 
densest peaks created by the turbulent compression. The 
mass in stars is still tiny, well under 1% of the gas mass, 
and the stars themselves are all quite small. Nonethe- 
less, the effects on the temperature are immediately ap- 
parent. Each star is surrounded by a clear region of gas 
at elevated temperatures. These regions are localized, 
so that the the bulk of the gas remains cold, and the 
heated regions around different stars are, for the most 
part, non-overlapping. 

It is not surprising that the format ion of stars has such 
a strong effect. As pointed out by Offner et al. (2009), 
the energy budget of a star-forming cloud is dommated 
by the gravitational potential energy released by star for- 
mation, even when those stars constitute a tiny fraction 
of the total mass. This continues to be true up until the 
point when massive stars with short Kelvin times begin 
to dominate the bolometric output of the stellar popu- 
lation. In our simulations, even though we do produce 
~ 20 Mq stars with significant internal luminosities to- 
ward the end of the simulations, accretion luminosity is 
the dominant energy source over most of the simulation 
time. 

This morphology of small regions of warm gas strung 
out along filaments continues to hold to some extent even 
at time t/tff = 0.6, when the stellar mass has increased 
to a few percent of the gas mass. We can still identify 
distinct heated regions associated with individual stars 
or small stellar groups, and the bulk of the mass remains 
near 10 K. In the last two time slices, however, as a 
larger and larger fraction of the cloud mass is converted 
into stars, this ceases to be true. Even the coldest gas 
anywhere in the cloud is now at temperatures notice- 
ably larger than the original background temperature, 
and the regions of very warm gas, T ^ 100 K, are begin- 
ning to overlap and merge. In the last time slice, the 
coldest gas anywhere in the computational domain is at 
~ 30 K, and much of the mass is concentrated in a few 
compact regions where the temperature is significantly 
higher. Rather tha n a few war m, dense regions around 
individual stars (cf. Offner et al. 2009) the bulk of the 
gas is now concentrated into a smaller number of more 
massive regions that are heated by the collective effects 
of large numbers of stars. 

3.2. Star Formation History and IMF 

Figure |4] shows the total mass of all stars as a function 
of time in the runs. Examining the figure shows that the 
total mass in stars is nearly identical in the two radia- 
tive runs, indicating that this aspect of the simulations is 
very well converged. Run ISO begins to form stars some- 
what earlier, and the mass in stars present at equal times 
is somewhat higher. However, this difference mostly ap- 
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Fig. 4. — Total mass in stars M*,tot, normalized to the initial 
cloud mass Mc, as a function of time, t, normalized to the mean- 
density free-fall time fff. We show results for run LR (thin red 
line), run HR (thick blue line), and run ISO (thin green line). 




Fig. 5. — Total number of stars (top) and number of stars 
normalized to the number present at the time when ^tot / Mc = 
0.5 (bottom). We show results for run LR (thin red line), run HR 
(thick blue line), and run ISO (thin green line). 

pears to be a time offset. The overall shape in Figure [4] 
is the same, indicating a generally similar star formation 
history. The time offset is likely a result of the faster 
collapse that occurs in the isothermal run, where cooling 
is assumed to be infinitely rapid and efficient, compared 
to the radiative run. 

Figure |5] shows the number of stars as a function of the 
total stellar mass in each simulation. The total number 
of stars is somewhat larger in run HR than in run LR, 
which is not surprising given the increased resolution. 
Observations indicate that the binary period distribution 
is extremely broad, coveri ng separations from only a few 
stehar radii to > 10"^ AU ( [Duquennoy fc MayorllOOTI ). It 
is therefore not surprising that some binaries that might 
be resolved into two separate stars in run HR instead 
appear as a single star in run LR - indeed, we would 
expect this result in essentially any simulation that did 
not resolve the radii of individual stars. Nonetheless, 
notice that, if we normalize to the number of stars present 
at equal times and fractions of mass accreted, then the 



difference between the two runs disappears. The number 
of stars present at any given time in run HR is roughly 
1.6 times the number present at the same time in run LR. 
Thus the trend in terms of when the stars are formed in 
the simulations is nearly identical in the two cases, and 
we can regard as well-resolved the distribution in time of 
when stars form. 

The trend of number of stars versus mass shown in 
Figure [5] is interesting. In the radiative runs, when 
M^^tot/Mc ^ 0.1, the number of stars increases roughly 
linearly with the total stellar mass, as we might ex- 
pect if the mass per star were constant. However, the 
rate at which new stars appears drops sharply once 
M^^tot/Mc > 0.2. Indeed, we see that 60 - 70% of ah 
stars have formed at a time when only ~ 10% of the 
cloud mass has been incorporated into stars. By the time 
20% of the cloud mass has gone into stars, nearly 90% of 
all the stars are in place. In effect, the fragmentation of 
the gas into new stars has completely shut down. Given 
that this effect occurs nearly identically in runs LR and 
HR, this cannot be a resolution effect. In contrast, run 
ISO shows very different behavior. The number of stars 
as a function of total stellar mass is almost the same as 
in run HR up to the point where ~ 15% of the mass 
has been incorporated into stars, but the two runs di- 
verge after that. New stars continue forming all the way 
through run ISO, at a rate that is only slightly less after 
M^^tot/Mc ^ 0.2 than it was earlier in the simulation. 
This strongly suggests that the shutdown in new star 
formation we observe in runs LR and HR is a radiative 



effect, a topic to which we will return in Section 3.3 



As one might expect, this shutoff of fragmentation into 
new stars in runs LR and HR even as the total stellar 
mass continues to increase produces a dramatic effect 
on the stellar mass distribution. Figures ^ and [7| show 
the cumulative and differential mass distributions of the 
stars formed in our simulations at the times when the 
total mass in stars is 10 — 50% of the initial cluster mass. 
All these plots show that the stellar mass distribution in 
the radiative runs moves continuously to higher masses 
as the simulation proceeds. This is because mass is ac- 
creting onto existing stars, which rise in mass, but very 
few new, lower-mass stars are forming. Note that, while 
the mean stellar masses are slightly different in runs LR 
and HR, the systematic drift of these mean to higher 
masses as the total stellar mass rises appears to about 
occur equally in both runs. In run ISO, on the other 
hand, there is much less evolution in the shape of the 
IMF. The fraction of mass in very small objects does de- 
crease slightly with time, but the IMF in run ISO peaks 
at ~ 1 Mq in every time slice. Quantitatively, we find 
that, from the point where M*^tot/^c ^ 0.15 and the 
star formation histories in runs ISO and HR begin to di- 
verge, up to the point when M^^tot/^c ^ 0.5 and run 
HR ends, the mass- weighted median stellar mas^in run 
ISO increases by only a quarter of a dex, while in run HR 
it increases by half a dex. Thus the behavior of run ISO 
is similar to that in previous simulations done with pre- 
scribed equations of stat^ (e.g. see Figure 1 of ,Bonnell| 



^ Defined as the mass m for whic h stars with masses m^c 
com prise half the total stellar mass 



< m 



' We cann ot directly compare to the earlier radiative simula- 
tions ot low mass clusters by Bate ( 2009 ) and Ottner et al. (2009 ), 
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Fig. 6. — Fraction /(< M*) of the total stellar mass in stars 
with masses smaller than M* , as a function of M* . The five panels 
show this cumulative mass function in the simulations at the times 
when the total stellar mass is 10%, 20%, 30%, 40%, and 50% of 
the initial cloud mass as indicated. We also show the times at 
which simulations LR and HR (but not run ISO) reach these stellar 
masses, which are identical to within a few percent in runs LR and 
HR. We show run LR (red line), run HR (blue line), run ISO 
(green line) and the results of creating a cluster of star s of equal 
total mass (gray), with each star randomly drawn from a|Chabrier| 
|2005) IMF following the procedure outlined in Appendix [a] tor 
the Uhabrier (2005 ) IMF, the three hues show the 10th, 50tlT; and 
90th percentile ot the random drawings, and the hatched region 
indicates the range between the 10th and 90th percentiles. 



et al. (2004), which a similar increase in median mass 



from 0.7 — 1.0 free-fah times in their simulation.^ 

For comparison, we have generated 10,000 clusters each 
of mas s 100, 200. 300, 4 00, and 500 M©, randomly drawn 
from a Chabrier (2005) IMFPlwith a minimum mass of 
0.05 Mq and a maximum 01150 Mq. We properly ac- 

because these produced fewer than 20 objects. Their IMFs are 
therefore much too sparsely sampled for it to be possible to make 
any meaningful statements about their time-dependence. 

^ There may also be difference between our simulations and those 
of Bonnell et al. due to differences in initial conditions (partly 
centrally condensed for us versus uniform density for them) and 
equation of state (isothermal for our run ISO, barotropic for them.) 

^ The argument has been advanced in the literature that the 
observed IMF in clusters with masses as small as a few hundred 
Mq is is truncated at high masses compared to Chabrier or similar 
IMF ( Weidner et al. (2010 j, but see Lamb et al. (2010), Calzetti 
|et al.pUlUj and 1^'umaga lli^et al.| ( |20lT[ ) tor observational counter- 
arguments). Here we are mostly interested in the peak of the IMF, 
not the high mass end where our simulations have too few stars to 
make statistically strong statements. We therefore proceed with 
the simplest assumption that there is no high mass truncation, 
since it makes no difference for our purposes. 



Fig. 7. — Same as Figure [6] except that we show the differen- 
tial mass distribution (i/(M* j/dlog M*, i.e. the value in each bin 
indicates the fraction of all stellar mass that lies in that bin. We 
normalized our distributions so that the sum of all bins multiplied 
by the bin width equals unity. As in Figure [g] we show run LR 
(red line), run HR (blue li ne), ru n ISO (green Ime), and clusters of 
equal mass drawn from a [Chabr ier (2005) IMF following the pro- 
cedure in Appendix (grayj. I'he histogram values we plot are 
the 50th percentile oiour experiments, and the error bars indicates 
the range from the 10th to the 90th percentile. 



count for finite sampling using the procedure described 
in Appendix |Aj As the plots show, the mass distribu- 
tion of stars formed in the radiative simulations drifts 
to systematically higher masses than the observed IMF 
once ~ 30 — 50% of the mass has been turned into 
stars. The disagreement is highly significant, and oc- 
curs at stellar masses that are extremely well-resolved in 
the simulations. For example, consider Figure [7| at the 
time when M*^tot/^c = 0.5. For run HR at that time, 
the mass in almost every bin from 1 — 10 Mq is above 
the 90th percentile of random drawings from a Chabrier 
IMF, while the mass in almost every bin below 1 Mq 
is below the 10th percentile of random drawings from a 
Chabrier IMF. Indeed, a Kolmogorov-Smirnov compari- 
son between the mass functions produced in the simula- 
tions and the Chabrier IMF shows that, with the excep- 
tion of the HR run at the point when tot/^c = 0.3, 
all the mass functions shown in Figures |6] and [7| are incon- 
sistent with having been drawn from the Chabrier IMF 
at confidence levels better than 1 part in 10^. 

In contrast, run ISO is consistent with the IMF at the 
low mass end at essentially all times. It is deficient in 
massive stars compared to a Chabrier IMF, an effect 
that has been observed before in simulations without ra- 
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diative transfer ( [Maschberger et aL 2010) and taken as 
evidence for the so-called 'IGIMI^' (integrated galactic 
initial mass function) effect". We obtain the same re- 
sult here, but find that it disappears in simulations that 
include radiation. 

One might be tempted to fix this problem simply by 
scaling all the stellar masses by some factor less than 
unity, to account for mass ejected by protostellar out- 
fiows, which we have not included. However, because 
the peak of the IMF is evolving with time in our simula- 
tions, a scaling factor that produces agreement between 
the simulated IMF and the observed one at one time 
would not produce agreement at earlier or later times. 
The central problem is not so much that the IMF in the 
simulation is too top-heavy, but that the median mass 
increases continuously with time. However, it does seem 
likely that protostellar outfiows can help solve the prob- 
lem by reducing the star formation rate and thus the 
luminosity, as we discuss further in Section 4.3 Such an 
effect cannot be captured by a simple rescaling of the 
masses, though. 

3.3. Gas Thermodynamics and Fragmentation 

The reason for the shutdown in fragmentation and the 
drift to systematically higher stellar masses with time in 
the radiative runs becomes clear if we consider how the 
gas density and temperature evolve with time. Figures [S] 
and [9] show the distribution of gas mass in the density - 
temperature plane as star formation proceeds in runs LR 
and HR. For comparison, we also overlay lines of constant 
Bonnor-Ebert mass, where 



Mbe = 1.18 



(15) 



and Cs = y^ksT/ ji is the isothermal sound speed. The 
Bonnor-Ebert mass is significant because objects with 
masses below Mbe can be supported against collapse by 
thermal pressure. We therefore expect that the lowest 
mass stars to form will tend to have masses comparable 
to the smallest values of Mbe found in the gas. Even 
if turbulence does create fragments with masses below 
Mbe, these will be stable against collapse as a result of 
their thermal pressure. Figure [To| summarizes this result 
by showing how the gas mass isdistributed with respect 
to Mbe at different times in the simulation. As the plot 
shows the runs are not completely converged at the low 
Mbe end, but the general trend that the mean Bonnor- 
Ebert mass systematically increases is clear in both runs. 

Figures |8]- [To] show that, immediately before any stars 
have formed, the great majority of the mass has a tem- 
perature within a few K of 10 K, the initial tempera- 
ture and the temperature imposed by the background 
radiation field. Consequently, the densest material in 
the cloud is in a density and temperature regime where 
Mbe ^ 0.01 M©, and objects of this mass are able to 
collapse. Nearly half the mass in the cloud lies in the 
region between Mbe =0.01 M© and 0.1 M©, so there 
is plenty of material available to make low mass stars. 
Once stars begin to form, however, their radiation raises 
the temperature significantly, pushing to higher Mbe- 
This increase is partly offset by an increase in the mean 
density, but the density does not increase quickly enough 
to compensate for the rapidly rising temperature - likely 



because the density rise occurs on a timescale associated 
with the mean density free-fall time, while the temper- 
ature rise is driven by stellar accretion occurring at the 
peaks of the density distribution, which operate on a 
much shorter timescale. As a result of this evolution, 
there is not much material that is dense and cold enough 
to make small stars. For example, if we look run at HR, 
we find that 20% of the gas mass has Mbe < 0.05 Mq 
just before the first star forms, and thus is able to make 
the smallest stars we consider. In contrast, the mass 
fraction able to create such small stars drops to 10% at 
the point when M*^tot/Mc = 0.2, and to only 2% when 
M*^tot/Mc = 0.5. Thus we see that the formation of new 
stars has stopped because it is no longer possible for small 
fragments to gravitationally collapse. By the end of the 
run, the smallest gravitationally-unstable fragments are 
approaching 1 Mq in size. 

The underlying physical reason for this effect, of 
course, is the radiation released by the already- formed 
stars. This in turn, is primarily driven by accretion lu- 
minosity, with a subdominant contribution from nuclear 
burning and Kelvin-Helmholtz contraction later in the 
simulation. 

4. DISCUSSION 
4.1. The Overheating Problem 

A systematic increase in the mean stellar mass induced 
by heating of the gas due to accretion luminosity is a 
new phenomenon in simulations of star cluster forma- 
tion. Radiative suppression of fragmentation has been 
reported in the literature before, but no previous sim- 
ulation has observed it to shift the typical stellar mass 
scale in regions as large as entire star clusters. We em- 
phasize that, even if we regard the absolute stellar mass 
peak we obtained as uncertain due to resolution effects, 
the trend of increasing mean mass is robust, and appears 
equally strong in both simulations. It has not been seen 
in previous work due to the limitations we outlined in 
Section [T] Most simulations of large-scale cluster for- 
mation with initial conditions similar to ours, which are 
typical of Galactic star formation, have not included ra- 
diative transfer. They have adopted a simple equation 
of state, which puts in by han d the result that the peak 
stellar mass is invariant (e.g. Bonnell et al. 2004). We 
effectively recover the same results in our run ISO: the 
median stellar mass does increase slightly with time, but 
the increase is significantly smaller than in the radiative 
runs, and is consistent with what has been found in ear- 
lier non-radiative simulations. 

Those simulations that have included radiation have 
either focused on refflons too small or wi th too few stars 
(e.g. Bate 2009; Offner et aI]|2009||2010D to produce the 
overlap of the heating regions around many stars we ob- 
serve, or have focused on single massive co res, where 
suppr ession of frag mentation is expect ed (e.g. Krumholz 
et al.||2007al 12010 iMyers et al.||201lD. Indeed, in these 



contexts, radiative suppression of fragmentation is neces- 
sary to obtain agreement between simulations and obser- 
vations. For single massive cores, suppressed fragmen- 
tation has tentatively bee n seen in high resolut i on in- 
ter ferometric obs ervations ( Bontemps et al.||2010 Long- 
more et al. 2011) . In the absence of radiation, the disks 



formed in simulations of low mass star formation tend to 
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Fig. 8. — Distribution of the gas mass in bins of density and temperature in run LR. The panels show the distributions at the time 
when the total stellar mass normalized to the initial cloud mass is M^^tot/Mc = 0.0, 0.1, 0.2, 0.3, 0.4, and 0.5, as indicated. For the time 
M*,tot/Mc = 0.0 we select the last time snapshot for which M*,tot = 0. Within a panel, the color of each pixel indicates the density of 
mass in that pixel, measured in Mq per dex in density per dex in temperature. The sashed black line labelled "sink creation" indicates 
the locus of density and temperature for which a computational cell exceeds the sink particle creation condition, equation (|14| evaluated 
with J = 1/4. Gas cells that fall to the right of this line create sink particles that incorporate some of their mass, which is why no cells fall 
to the right of the line. For comparison, the dotted lines indicate the loci where the Bonnor-Ebert mass Mbe = 0.01, 0.1, and 1 Mq, as 
indicated. 

has been incorporated into stars. Figures [6j and [7| show 
that the shift of the simulation IMF to higher masses 
than the Chabrier (2005) IMF is also largely complete 
by this point. Since this is about the minimum star 
formation efficiency req uired to have any possibility of 



undergo excessive fragmentation, leading to an overpro- 
duction of brown dwarfs relative to stars and to various 
other confficts with observation (Luhman et al. 2007). 
Inclu ding ra diation fixes thi s problem (^Bate,2009^; Ottner, 
et al.||2009[|2010l) . Indeed, [Bate! ( |2009p argues that the 
observed peak or the IMF can be explained as arising 
from the mass scale at which radiative feedback halts 
fragmentation. While this argument is plausible, it re- 
lies on the assumption that we can consider the bubble 
of radiatively-warmed gas around each star to be iso- 
lated from other stars, amidst a background of cool gas. 
This assumption hol ds in the l ow- mass, low-dens i ty re- 
gions simulated by |Bate| (|2009) and [Offner et aD ( |2009 , 
2010), where regions of heating are ^ 0.05 pc m size, 
much smaller than the interstellar separation. It clearly 
does not hold in our simulation, both because our stars 
are closer together than in a low mass star-forming re- 
gion, and because our heating regions are larger due to 
the higher accretion rates produced by the higher gas 
densities. This suggests that the critical problem in our 
simulation is that the regions of warm gas around indi- 
vidual stars begin to overlap. As a result, all the gas in 
the cluster is heated, rather than simply discrete regions. 

One might hope to avoid this problem by halting star 
formation early on, before enough mass goes into stars 
to allow the heated regions to overlap. However, such 
a solution seems to require improbable fine tuning. Ex- 
amining figures [T] and |2j we see that the overlap of hot 
regions is well underway by the time 30-40% of the mass 



making bound clusters (Kroupa et al. 2001; Fall et al. 
2010), the fact that at least some star lorrnation does 



result in bound clusters suggests that the star formation 
efficiency cannot be vastly lower than this value most of 
the time. 

4.2. Understanding the Problem 

We can estimate the dividing line between the two 
cases of heating in discrete regions around single stars 
and heating in the bulk of the protocluster gas using the 
analytic radiative tran sfer approximation of Chakrabarti] 
fc McKee|(|2005[ 120081), coup l ed to t he formalism devel- 
oped by |Krumholz <^ lV[cKee| ( |2008[ ). [Chakrabarti fc Mc- 
Kee consider a spherical cloud of dusty gas with radius 
i?, mass M, and a density profile p oc r~^^, surround- 
ing a point source of radiation of luminosity L, with 
dust whose specific opacity depends on wavelength as 

= hZo{Xo/X)^. In such a cloud, they show that the 
temperature profile approximately follows 



T = Tch 



R 



ch 



-kT 



(16) 



where r is the distance from the cloud center, R^h and 
Tch are the characteristic radius and temperature of the 
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Fig. 9. — Same as Figure [8] but for run HR. Note that the sink creation locus is shifted to higher densities by factor of 4 compared to 
run LR. 




Fig. 10. — Fraction f{M < Mbe) of the gas mass for which 
the Bonnor-Ebert mass a t th e local density and temperature, as 
computed from equation JTsJ, is less than Mbe- We show these 
mass distributions for run LR (red) and run HR (blue), at the 
times when the total mass is M^^tot/Mc = 0.0 (dashed, just before 
the first star forms), 0.2 (dot-dashed), and 0.5 (solid). 



dust photosphere formed within the cloud, and kr is a 
powerlaw index to be approximated by a numerical fit. 
For convenience we define S = M/irR^, rj = L/M (mea- 
sured in cgs units, not Solar units), a = l/[2/3+4(/cp — 1)], 

R R/Rch, and To hc/kBXo- For Milky Wa y dust, 
0^2 and ^ 0.54 c m^ g"^ at Aq = 100 fim (|Wein- 



gartner fc Draine||2001| ), but the results depend on these 
parameters very weakly. With these definitions, the char- 
acteristic radius and temperature and the powerlaw in- 
dex are given by 



-^ch 

~R^ 



^cr SB L 



3 - kp)f^o 
4{k, - l)T^ 



(17) 



)0.1 



4(fcp - 
(3 - kp)Ko 



0.48fc°°^ 



R 



(18) 

(19) 
(20) 



The latter two expressions are approximations based on 
fits to numerical solutions of the radiative transfer equa- 
tion, and reproduce the numerical results with high ac- 
curacy. 

A rough condition for the heating regions around indi- 
vidual protostars to merge and heat the bulk of the gas 
is that the combined luminosity L of all the protostars, 
which we approximate as being near the cloud center, be 
high enough so that the temperature T at the edge of 
the cloud, r = i?, be higher than the background tem- 
perature T5 « 10 K to which gas settles when it is not 
heated by a nearby star. Thus, to avoid overheating we 
require that the luminosity to mass ratio 77 be smaller 
than the value r^crit for which T{R) = T5. For a given 
cloud mass M and surface density S, it is straightforward 
to use equations (16) - (20) to numerically determine the 
value 77crit for whicn the condition T{R) = T5 is satisfied. 
In what follows we do so for a background temperature 
T5 = 10 K and density profile /Cp = 3/2, roug hly what is 
seen in massive star-forming clumps (e.g. Mueller et al.| 
2002), but the result we obtain is not very sensitive to 
this choice. 

The luminosity is in turn related to the star forma- 
tion rate in the simulations. Krumholz & McKee (2008) 
show that accretion onto low mass stars yields an en- 
ergy release per unit mass accreted i/j ^ 2.1 x 10^^ erg 
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Fig. 11. — Luminosity to mass ratio 77 versus dimensionless 
star formation rate Cff. In each panel, dashed horizontal lines 
show the critical light to mass ratio ?7crit at which heating re- 
gions around individual protostars merge and fragmentation is sup- 
pressed. Solid lines show the value of ?7acc for an accretion-powered 
cluster-forming gas cloud with the indicated value of Cff. In the 
upper panel we show clouds with E = 0.5, 1, and 2 g cm"-*^ (red, 
green, and blue), all of mass M = 10^ Mq. In the lower panel we 
show clouds of fixed E = 1 g cm"-^, with masses M = 10^, 10^, 
10^, and 10^ Mq (red, green, blue, and purple). 

g~^. This number is a result of stellar structure consid- 
erations, which fix the characteristic radii of protostars. 
Thus the light to mass ratio T^acc in a cloud of mass M 
powered by accretion luminosity from stars forming at a 
rate is 



M 



eff 



(21) 



where tf[ = ^/ttWJsGM is the mean-density free-fall 
time of the cloud and e g is the dimensionless star fo rma- 
tion rate introduced by [Krumholz 



(2005) 



McKee| 

In Figure [11] we plot r/crit and 77 lor clouds of varying 
mass M and surface density S as a function of e^. Our 
values of M and H are chosen to span the range of typ- 
ical cluster-forming gas clumps in the Milky Way. The 
value of T^crit of course depends on S alone, while that 
of T^acc is proportional to eff. The plot shows that, for 
any plausible cloud mass and surface density, if > 0.1 
then 77 > T^crit- This plot explains why, in our simulation, 
the stellar heating regions overlap. In our simulation, 
~ 50% of the gas is in stars when t/tff 1, so we have 
eff ^ 0.5. This puts us in the regime in which heating 
zones overlap, and fragmenta tion is suppressed. In con - 
trast, the simulations of Bate] ( [20091 ), [Off ner et al.| ( |2009 D, 
and Pete rs et al.| ( |2010p have substantially lower surlace 
densities, S < 0.1 g cm~^, putting them in the regime 
where heating zones do not overlap and fragmentation 
is unlikely to be suppressed even with quite high e^, ex- 
ce pt w ithin the disks around each star. The simulations 
of |Offner et al. , since they include driven turbulence, also 
have a lower value of . 

Note t hat this argument is c onsistent with the one 
made by Elmegreen et al. ( 2008 ) for why the Jeans mass 
should vary little in regions where the dust and gas tem- 
peratures are well-coupled, like those we study. The crux 
of their argument is that increases in the gas density 



lower the Jeans mass, but also produce a higher star for- 
mation rate, which in turn raises the dust temperature 
and the Jeans mass. The two effects nearly offset one 
another. However, this offset only occurs if the star for- 
mation rate and the density are related by a volumetric 
Schmidt (1959^ .1963 ) law with eff ~ 0.01 (see their equa- 
tion 6). It eff rises with time, as it does in our simulation, 
then the Jeans mass will not remain independent of den- 
sity. 

4.3. Possible Solutions to the Problem 

Unders tanding the problem also suggests an immediate 
solution. Kru mholz fc Tan] ( 2007 ) compile observations 
from the literature and tina th at, even i n dense, cluster- 
forming gas clumps, eff ^ 0.01. jEvans et al. (2009) obtain 



similar values of eff in cluster-forming regions observed 
by the c2d Survey, although c2d targets clusters consider- 
ably more diffuse and lower in mass than the one we have 
simulated here. Figure 11 shows that such clouds are not 
in the regime where heating zones overlap and fragmen- 
tation is suppressed, unless their masses are quite low, 
M < 10^ Mq. This explains why fragmentation is not 
suppressed in real clouds. 

However, this result has important implications for 
simulations of star cluster formation. It implies that, 
once radiation physics is included in the simulations, one 
cannot expect to obtain the correct IMF without also ob- 
taining the correct star formation rate, or at least a star 
formation rate that is roughly correct. In simulations 
that do not include radiation feedback, no such care is 
needed. Fortunately, obtaining the correct star forma- 
tion rate in simulations is not terribly difficult. Sim- 
ulations where the turbulence is driven artificially can 
reproduce the observed value ^ 0.01. Even better, 
simulations that include stellar wind feedback naturally 
produce realistic, l ow values of eff without any artificial 



driving (e.g. Li fc Nakamura ^2006VNakamur a fc Lil2007 
Wang et al. 2010l), and preliminary evidence indicates 
that this does reduce accretion luminosities to the point 
where fragmentation is suppressed far less than we have 
found (Hansen et al., 2011, in preparation). It is not 
clear if this effect is scalable to all cluster masses ( Fall] 
et al. |2010 ), but it does suggest a way toward simula- 
tions of cluster formation that simultaneously obtain the 
correct star formation rate and the correct IMF. 

It is thus possible that the problem might be solved the 
the inclusion of other physics that our simulation omits, 
such as outflows and photoionization. These mechanisms 
might be able to generate regions of high enough density 
that their Bonnor-Ebert masses will be low even in the 
presence of overlapping regions of radiative heating. 

4.4. Implications for Competitive Accretion versus Core 

Accretion 

It is also interesting to consider the implications of our 
result for the competitive accretion versus core accre- 
tion models for the formation of star clusters and origin 
of the IMF. Roughly speaking, the competitive accre- 
tion model is that collapses that produce star clusters 
are global in nature, so all stars accrete from the same 
mass reservoir, and the stellar mass distribution is de- 
termined by a competition between formation of new, 
small fragments (which pushes the mean mass to lower 
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values) and growth of existing fragments by Bondi-Hoyle 
accr etion (which pushes the mean mass to higher values) 



(e.g. i Bonneh et al.|[2QQT| [Bate fc Bonnell||2QQ5[ iBonneh 
|& Bate 2006). in contrast, in the core accretion model, 



collapses that produce individual star systems are local 
rather than global, so that different protostars are for the 
most part not accreting from the same mass reservoir. In 
this case, the mass distribution of the stars is set by the 
mass distrib ution of the regions of localized collapse, the 



cores 



2003 



Padoan fc Nordlund||200 2; McKee fc Tan 
Fado an et al.||2007| |Alves et al.| ,2007; HennelDelle 
labrier 2008; 120091). Intermediate models are also 



possible, in which low mass stars form via local collapse, 
but either massive stars or the cor es from which they 
grow form via a global collapse (e.g. Peretto et al.||2006 
Wang et a l. 2010). 
[KrumEolz^et al. (2005|) point out, andlBonnell & Bate 



(|2006|) and |Ottner et al.f (|2008a|) confirm, that whic. 
mode of star formation takes place depends on the level 
of turbulence and on Cff. If the turbulence is sub-virial, 
or becomes sub-virial through decay that is not offset by 
internal feedback or external driving, then eff becomes 
large and competitive accretion is the dominant star for- 
mation mode; core accretion prevails if the turbulence 
remains at virial levels and is small. In our simu- 
lation we do not include either artificial driving or any 
physical feedback mechanisms capable of driving the tur- 
bulence (e.g. protostellar outfiows or H ll regions), so 
our simulation produces large and we obtain a com- 
petitive accretion-like mode of star formation. However, 
crucially, we have shown that such a mode of star forma- 
tion cannot produce the correct IMF due to the radia- 
tive suppression problem we have identified. The con- 
stant production of new, low-mass stars on which com- 
petitive accretion relies to keep accretion onto existing 
stars from pushing the IMF to ever-increasing masses 
does not happen once radiative feedback is included, at 
least in the minimal case where hydrodynamics, gravity, 
and radiative feedback are the only physical ingredients. 
It is conceivable that some mechanism we have omitted 
might still enable the production of low mass stars even 
in clusters with high eff (e.g. fragmentation induced by 
expanding H ll shells), but in this case that mechanism 
would be responsible for controlling the peak of the IMF. 
Our results therefore suggest the minimal competitive ac- 
cretion model is not compatible with the observed IMF. 

One might try to alleviate this problem by choosing 
significantly less dense initial conditions and retaining 
the high required for competitive accretion, i.e. by 
selecting a lower S in Figure [TT] However, this solution 
faces a severe problem: the initial conditions we have 
selected are typical of the observed gaseous properties of 
clouds wher e massive star formation occurs (e.g . Shirley 



et al. 



2003 



Faiindez et al, 



2004 Fontani et al. 2005). 



JSurface densities are even larger m globular clusters, yet 



these show the same IMF peak as the field De Marchi 
et al. ( 2000[ [2010). If the minimal competitive accretion 



model can only reproduce the observed IMF from initial 
conditions far less dense than we have simulated, then its 
applicability is limited to low-density regions like Taurus, 
which generally do not contain any massive stars. 

4.5. Implications for Fragmentation- Induced Starvation 




Fig. 12. — Stellar mass as a function of time for a sample of 
individual stars in runs LR (top), HR (middle), and ISO (bottom). 
The stars shown are the five most massive at the final time in the 
simulations, plus 10 other stars evenly distributed in mass. 



10"^ 



10"^ 




Fig. 13. — Mass accretion rate for the most massive (solid) and 
second most massive (dashed) stars at the final time in runs LR 
(red), HR (blue), and ISO (green). To minimize confusion, the 
accretion rates have been smoothed over a timescale of O.OStff. 

It is also interesting to examine how individual stars, 
and particularly the most massive stars, grow in mass. 
We show this in Figures 12 and[T3j which show the mass 
versus time and the mass accretion rate versus time for a 
sample of stars in each run. In run LR the most massive 
star we form is 20.0 M©, in run HR the most massive 
star is 16.2 Mq, and in run ISO it is 10.3 Mq. In runs 
LR and HR, the most massive stars are continuing to 
grow rapidly at the end of the simulation, with accretion 
rates that are generally fiat or increasing with time. The 
most massive stars are also growing in run ISO, but more 
slowly and with accretion rates that are either constant 
or declining with time. 

These results have potential implications for t he idea 
of fragmenta tion-induced starvation proposed by Peters 
et al. (2010'). In their simulations (which have a resolu- 
tion comparable to that of our run LR) , the most massive 
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stars stop growing after a certain point because the ac- 
cretion flow that is feeding them fragments to produce 
smah stars rather than being accreted by the massive 
star. Thus massive stars exhibit accretion rates that fall 
with time. We do so something roughly consistent with 
this behavior in run ISO, but not in our radiative runs. 
This is likely an effect of radiative suppression of frag- 
mentation. 



Peters et al. (2010) also include radiative transfer in 
their simulations, but they do not find strong suppression 
of fragmentation. This is probably because their simu- 
lated cloud has a much lower column density (S ~ 0.03 
g cm~^) than either our simulated clouds or than typical 
regions o f massive star for mation in the Galaxy (S ~ 1 



g cm ). Krumholz et al. (2010) show that the amount 



by which radiation suppresses fragmentation is highly 
sensitive to the column density, and predict essentially 
no su ppression at the column density used by Peters^ 



et al. The physical reason for this is that a cloud with 
S = 0.03 g cm~^ is optically thin even in the near- 
infrared, so starlight that is absorbed by dust grains 
promptly escapes, and most gas is not heated by the radi- 
ation. It is therefore not surprising that Pe ters et al.^ see 
fragmentation- induced starvation and we do not. 

We emphasize, however, that the absence of 
fragmentation-induced starvation in our radiative runs 
does not mean that fragmentation-induced starvation 
does not occur under typical Galactic star-forming con- 
ditions. We have just argued that fragmentation is sup- 
pressed too strongly in our simulations because star for- 
mation is too rapid. Indeed, simulations indicate that 
outflows allow more fragmentation to occur even in sin- 
gle massive cores than in comparable simulations with- 
out outflows (Cunningha m et al.||2011| . However, our 
results suggest that, before tragmentation-induced star- 
vation can be considered an important mechanism in 
regulating massive star formation, it will be necessary 
to simulate the formation of a star cluster using typical 
Galactic conditions, like we do, and to include mecha- 
nisms that produce realistically low star formation rates. 

5. SUMMARY 

We report simulations of the formation of a mas- 
sive star cluster comparable in size to the Orion Neb- 
ula Cluster. Our simulations use adaptive mesh refine- 
ment to obtain high resolution, and include radiation- 
hydrodynamics coupled to a realistic treatment of stellar 
radiative feedback. These are the first simulations re- 
ported in the literature that include radiation feedback 
in the context of the typical region of Galactic star clus- 
te r formation, as opposed to focusing o n single low-mass 

f!ommerco n et alj|2010|) or hig h- mass (Krumholz et al. 
07a, 2010; Myers et al. 2011) cores, or on low-mass or 
v-density regions like laurus ( Bate|[2QQ9 Offner et al. 
20101 [Peters et al.|[2QTQl ). 

Our simulations return a surprising result. At early 
times in the simulations, accreting stars produce bubbles 
of warm, radiatively-heated gas around themselves, and 
within these bubbles fragmentation is suppressed by the 
increased Bonnor-Ebert mass. However, we find that. 



once ~ 10 — 20% of the gas in the protocluster has been 
converted to stars, these bubbles of warm gas begin to 
overlap and merge. Rather than resembling a few warm 
islands surrounded by a sea of cold gas, we instead have 
a cloud where all the gas is warmed by the collective 
luminosity of all the accreting stars. 

Once the simulation reaches this state, radiation feed- 
back raises the temperature and the Bonnor-Ebert mass 
throughout the remaining gas enough to essentially halt 
the formation of any further stars. Mass continues to be 
converted from gas to stars, but this is almost entirely 
through accretion onto existing stars rather than forma- 
tion of new ones. As a result, when radiation is included, 
the stellar mass distribution in a globally-collapsing star 
cluster such as the one we simulate is not nearly con- 
stant or very slightly increasing with time, as has been 
reported in earlier, non-radiative simulations, and as we 
find here in a control run that does not include radiation. 
Instead, the stellar mass distribution shifts strongly to 
systematically higher masses as star formation proceeds, 
eventually becoming too top-heavy compared to the ob- 
served IMF. While the absolute mass scale remains un- 
certain in our simulations due to our inability to resolve 
tight binaries, the result that the IMF is non-constant 
and increasing with time is robust against changes in 
resolution. This implies that, unless there is also some 
mechanism to ensure that star formation in every proto- 
cluster stops when the IMF peak is in the same place, it 
is not possible to produce the invariant IMF peak that 
we observe via the global collapse scenario we have sim- 
ulated. 

We argue that the underlying reason that this problem 
occurs is that, in the absence of either external turbulent 
driving or any sort of internal mechanical feedback to 
slow star formation down, stars in our simulation form 
too quickly. Since accretion luminosity produced as gas 
falls onto stars is what ultimately drives the temperature 
increase in our simulations that shuts off fragmentation 
and leads to a top-heavy IMF, the problem is likely to be 
alleviated in simulations that include enough physics to 
obtain a low star formation rate similar to that observed 
in real star clusters. We are in the process of conducting 
such simulations now, and will report on the results in 
future publications. 
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APPENDIX 

A. GENERATING COMPARISON IMF SAMPLES 

The statistical samples for the IMFs shown in Figures |6] - [7| consist of stellar populations drawn from a |Chabrier| 
([2005") IMF, subject to the constraint that the total mass of the population have specified value. We create each cluster 
by the following procedure. First, we draw stars from the C habrie r IMF, 

dn _ .rj exp(-[ln{M*/M0}-lnO.2]V2a2), < Mq .... 

d\nM, \ exp[-(lnO.2)V2a2](M*/M0)-l•3^ > Mq ' ^"^'^ 

where AT is a normalization constant and a = 0.55 In 10. We truncate this mass function at 0.05 Mq on the lower end 
(to match our minimum stellar mass in the simulation) and at 120 Mq on the upper end. We continue to draw stars 
so as long as the total mass of stars is smaller than the specified target mass. If we draw a star of a mass such that 
adding it to our population causes the total mass to exceed the target mass by more than 0.1 Mq, we reject it and 
draw another. We continue drawing until the total mass of stars is within 0.1 Mq of the target mass. 

Once we have a set of stars, we form the cumulative and differential distributions. We repeat this procedure 10,000 
times each for clusters of total mass (from 100 — 500 Mq. To produce the values shown in Figures |6] and l7| at each 
mass on the x-axis, we sort the values of the 10,000 cumulative or differential distributions at that value of M*. 
The 10th, 50th, and 90th percentiles shown are the values at that mass point or mass bin are the 1,000th, 5,000th, 
and 9,000th vales in the sorted lists. 
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